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This paper is devoted to the basic question of what factors determine the strain field in a quasi 
static granular flow. It is shown that using stress - strain rate relations is not the proper way to 
understand quasi static rheology. An alternative approach is discussed where the local deformation 
is regarded as the cause of deformation in the vicinity. We suggest a continuum model where the 
local shear strain is proportional to its Laplacian and the proportionality factor is determined by 
the local stress. The predicted behavior is tested in a three dimensional shear cell by means of 
computer simulations. The simplicity of our setup makes it ideal to demonstrate and examine the 
fundamental open questions of collective granular flows. The observed shear profile is interpreted 
in the framework of the suggested model. 



I. INTRODUCTION 

When a granular material is deformed slowly in a sim- 
ple plane shear cell [jHH a critical state is reached after 
a transient period. This critical state is unique for the 
given granular material and is characterized e.g. by a 
well defined critical volume fraction $ C rit and a critical 
coefficient of effective friction /x cr ;t. A nice property of 
the critical state is that it is devoid of memory effects 
in the sense that the preparation history does not play 
any role [U, H| ■ The initial state of the contact network is 
destroyed and a self-organized inner fabric is maintained 
during the shear deformation. 

The material constants $ C rit and /x cr ;t are valid for 
quasi static shear of low inertia numbers I [l|, i.e. when 
the shear rate is low enough and the pressure is large 
enough. In this paper we concentrate on the quasi static 
regime rather than on dense inertia flows. The dilute 
gas-like regime is out of scope of the present study. 

The idea in the plane shear experiment is that the ma- 
terial is tested under uniform stress and strain conditions 
in order to deduce constitutive relations between local 
strain, strain rate and stress. Then with this knowledge 
one can try to understand the rheology in more compli- 
cated situations that often come up in applications. It is 
expected that a small homogeneous piece of bulk material 
of a large system behaves according to the plane shear 
experiment. One of the most interesting things about 
quasi static deformations is that the above concept does 
not work. 

The meaning of the effective friction coefficient is the 
ratio of the shear stress measured in the shear direction 
divided by the perpendicular normal stress. What we 
learn from plane shear tests performed both in exper- 
iments and in computer simulations is that persistent 
shear flow needs a stress ratio at least /U cr it. Once the 
imposed stress ratio drops below this threshold the flow 
stops and the material behaves as a solid body. 
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On the other hand if a bit more complicated setup is 
considered, e.g. a Couette cell the above rule "solid- 
below-the-thrcshold" does not hold any more. The prob- 
lem is more pronounced for split bottom shear cells [I|-0] 
where wider shear zones and smoother flow profiles arise 
than in Couette flows. In these systems stationary shear 
flow is observed also in regions where the local stress ra- 
tio is far below the threshold /x C rit [1, An extreme 
example is shown in (Io| where sand behaves as a zero 
yield stress fluid, i.e. even a tiny shear stress is enough 
to cause deformation. 

We conclude that such a collective flow can not be 
interpreted locally as a simple 1 plane shear. Clearly, the 1 
local strain rate is not determined by the local stress ten- 
sor alone. What is then the underlying mechanism that 
determines the stress and strain fields? How does the 
material know locally how to deform? This is the basic 
mystery of quasi static shear flow that is in the focus of 
the present paper. To answer these questions is of funda- 
mental importance in order to be able to predict rheology 
in various granular systems. Even a very simple geome- 
try such as the straight split bottom cell 0, Q represents 
a serious challenge for theories. We are not aware of any 
model in the granular literature that is able to provide 
the correct flow profile for this case. Let us refer to the 
above questions as the problem of collective rheology. 

A popular approach to describe the flow is to assume 
the effective friction coefficient to depend on the local 
inertia parameter fi e g(I) and indeed it does a very good 
job for inertia flows [H.Tll|. This technique fails, however, 
for quasi static flows. It can not describe the non-trivial 
smooth flow that emerges in the collective rheology and 
often predicts infinitely narrow shear bands instead 0, 

To study the problem of collective rheology it is a good 
start to deal with stationary flows for the sake of simplic- 
ity. Thus Couette or split bottom tests seem to be good 
choices. However, these systems are much more com- 
plicated than the plane shear experiment, they involve 
too many degrees of freedom and unnecessary difficul- 
ties. Just to mention one, the material layers that slide 
next to each other during the flow are curved in one way 
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FIG. 1: Schematic illustration of the shear cell used in the 
computer simulations. The shear cell is three dimensional 
where Lees- Edwards boundary condition is applied in y direc- 
tion while in x and z directions normal periodic boundaries 
are used. 



or another, which affects the stress field in a non-trivial 
way [f| and makes the treatment harder. To shed more 
light on the mystery it is better to start with the simplest 
possible system where the problem of collective rheology 
arises. In this article we suggest and analyze such a test 
system, which is only one step away from the homoge- 
neous plane shear. 

The shear flows investigated here are achieved inside 
the computer by means of discrete element simulations. 
We deal with three dimensional flows of frictional and 
spherical grains. The paper is structured as follows. 
First we examine the important reference situation of the 
quasi static critical state produced in homogeneous plane 
shear. Then we modify the system to achieve non-trivial 
collective flow and present the results of the simulations. 
In the last part we suggest a model that can account for 
some of the features of this collective behavior. 



II. HOMOGENEOUS PLANE SHEAR 
A. The influence of the inertia parameter 

The computer simulations are carried out by the al- 
gorithm of contact dynamics (CD) [l3l - [l6j ]. The dynam- 
ics, as usual, is governed by Newton's equation of mo- 
tion applied for each grain. The interaction of the grains 
arc handled by means of constraint forces at the contact 
points which are determined based on the assumption 
that the grains are rigid (undeformable) . 

Our shear cell uses Lees-Edwards boundary condition 
in order to attain homogeneous plane shear (see Fig. [I}. 
That is the system is periodically connected in y di- 
rection in a special way: There is a velocity difference 
Wshear between the upper and lower boundary which pro- 
vides the shearing. In x and z directions normal periodic 
boundaries are used. Thus confining walls are completely 
avoided and shearing is achieved by the special periodic 
boundaries alone. This is advantageous because walls 



would introduce undesired boundary effects (e.g. layer- 
ing) that significantly can alter the behavior especially of 
small systems used in computer simulations. This way 
the system is a priory translational invariant in x, y and 
z directions, there are no special points, everywhere is 
"inside the bulk" . The pressure conditions are controlled 
by Andersen boundary conditions 17j, i.e. by properly 
varying the volume. Details how this is implemented in 
CD can be found in fl5j . In all the simulations gravity is 
set to zero. 

Thus our system is sheared in the xy plane in x di- 
rection and the two following parameters are kept con- 
stant during the simulation: the global shear rate 7 = 
w s hear/(2Ly), where L y is the system size in y direction, 
and the yy component of the stress tensor (cr vy ). All sim- 
ulations presented in this paper are performed with the 
help of the above shear cell. 

In the first set of simulations we examine one system 
in various conditions. This test system consists of 2500 
spherical grains with radii distributed uniformly between 
the unit length 1.0 and 1.3. The friction coefficient /-t CO nt 
of the grain-grain contacts is set to 0.2 (for both static 
and dynamic friction). We are interested in the steady 
state properties therefore relatively large shear displace- 
ment is simulated for each run. The value of the total 
shear strain 7 = AZ s h e ar/ (2£j,) (given by half of the total 
shear displacement in x direction divided by the system 
size in y direction) is typically of the order of 100. Even 
for the smallest inertia number where the simulation is 
very computation consuming 7 is chosen larger than 20. 

We examine whether the inertia parameter / [l[ serves 
as a good control parameter that determines the state of 
the system. The inertia parameter is defined as 

I = id\fpjp , (1) 

where d, p and p are the typical grain size, the mass den- 
sity of grains and the pressure, respectively. The inertia 
parameter reflects the scale of the inertia forces with re- 
spect to the force scale generated by the pressure. Low 
inertia parameter means that the inertia forces that cause 
accelerations of the grains are much smaller than the typ- 
ical contact forces between the grains. The definition of 
p is ambiguous because the stress tensor is not spheri- 
cal. Here we use a yy in place of p, however, this choice 
is unimportant in the scope of the present study. Using 
any of the normal stress components or Tr(cr)/3 would 
lead to negligible change in the value of /. 

We focus basically on two indicators of the state of the 
material, the density and the resistance against shear. 
More precisely, these are the solid fraction $ which gives 
the total solid volume of the grains divided by the volume 
of the system and the effective friction /j, e g provided by 
the ratio of two elements of the stress tensor a xy /a yy . Wc 
measure the time averaged $ and fj, e s in the stationary 
flow and plot these data as the function of the inertia 
parameter. 

In Fig. [2] the results for the solid fraction are shown. 
We scanned through a wide region of / by systematically 
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FIG. 2: The solid fraction $ is plotted against the inertia 
parameter /. The circles connected by the solid line show a 
set of simulations which differ only in the driving shear rate 
7, while squares stand for simulations where other parameters 
are also varied in a random manner such as the mass density of 
grains, the pressure and the technical parameter At which is 
the time step of the algorithm. The inset shows the difference 
from the critical value (t&crit — <&(/)) against I in log-log scale. 
The straight line indicates slope 1. 



changing the driving shear rate 7, while keeping all other 
input parameters fixed. Other simulations were also per- 
formed where we tested several random combinations of 
the input parameters: The pressure and mass density of 
the grains were changed by factors over 1000, the time 
step used by the simulation code was also varied by a 
factor 50. It can be seen that putting all theses data in 
one plot make them collapse on a single curve. The in- 
fluence of other physical and technical parameters (no t 
shown here) were also tested: The size of the system [2(| , 
the inertia used by the Andersen ty pe o f pressure control 
and the number of force-iterations [141 ] used by the con- 
tact dynamics algorithm are all parameters that within 
a wide range have no effect on the value of $. 

Thus we conclude that in case of our virtual material 
the inertia number / is indeed the key control parame- 
ter and it determines the solid fraction uniquely in ho- 
mogeneous plane shear. Exactly the same holds for the 
effective friction /x e ff and other ratios of the elements of 
the stress tensor which are also unique functions of / as 
shown in Fig. [3] and Fig. SJ 

In Figs. O EH and 21 where the / axis is in logarithmic 
scale, all the plotted quantities reach a plateau for small 
inertia numbers. This plateau region can be regarded ap- 
proximately as the quasi static region. Its upper bound, 
i.e. the largest quasi static value of /, depends on the 
required accuracy. The ideal quasi static flow is, in fact, 
a mathematical limit case of I — >• 0. The homogeneous 
and stationary plane shear in the quasi static limit can 
be characterized by the well defined critical values: 



FIG. 3: The effective friction /x e fi is plotted against the inertia 
parameter /. The symbols (circles and squares) have the same 
meaning as in Fig. [2] The inset shows the deviation from the 
critical value (/i e ff(7) — ^crit) against I in log-log scale. The 
straight line indicates slope 1. 
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FIG. 4: The ratios of normal stresses are shown as the func- 
tions of the inertia parameter I. Open and closed symbols 
stand for a xx /a yy and a zz /a yy , respectively. Circles and 
squares have the same meaning as in Fig. [2] 



In our simulations we find that the deviations of $ and 
jj, e ff from the corresponding critical values are propor- 
tional to / for small inertia numbers (insets in Figs. [21 
[3]). 

An interesting feature of the normal stresses is shown 
in Fig. |U Regarding the xy shear plane the values of nor- 
mal stresses a xx and a yy are very close to each other pro- 
vided the system is sheared quasi statically. In perpen- 
dicular direction, however, the normal stress a zz turns 
out to be lower about 10 percent. This weak stress re- 
sponse in side direction might be an important factor 
in real three dimensional flows and lead to a non-trivial 
stress distribution throughout the system. 

In this section we examined how the quasi static crit- 
ical state (in short critical state) can be approximated 
by low inertia numbers. The simulations presented in 
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FIG. 5: The solid fraction <£> is plotted against the friction co- 
efficient of the grain-grain contacts. The data shown by circles 
and squares are obtained in polydisperse and monodisperse 
systems, respectively. 
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FIG. 6: The effective friction fj, e g is plotted against the fric- 
tion coefficient of the grain-grain contacts. The data shown by 
circles and squares are obtained in polydisperse and monodis- 
perse systems, respectively. 



the rest of the paper remain as near as possible to this 
ideal quasi static limit. This effort is bounded by the 
available computational capacity because lower inertia 
numbers require more calculations. 



B. The effect of contact friction and size 
distribution on the critical state 

In this section we study how the critical state differs for 
different materials. The materials used in our simulations 
differ cither in the coefficient of contact friction /i or in 
the size distribution of the spherical grains. We use either 
a polydisperse system similarly to the previous section, 
i.e. with grain radii chosen uniformly between 1.0 and 
1.3, or a monodisperse system with grain radius 1.0. The 
inertia parameters are set to 5.2 * 10 -3 and 4.5 * 10~ 3 for 
polydisperse and monodisperse simulations, respectively. 
The systems consist of 5000 grains each and are tested 
with various values of \x. The properties of the critical 
state are measured and averaged over a typical total shear 
strain 7 « 30. 

In Fig. 03 [5] and [7J we show the behavior of the quan- 
tities 4>, /i e s and the normal stress ratios. Because the 
applied inertia parameters are relatively low, these mea- 
sured values can be regarded as approximate critical val- 
ues. Actually, the deviations between the data shown 
here and the exact critical values can be estimated based 
on the results of the previous section. For example the 
exact curve /i cr it (fi) is expected to be slightly lower than 
the values of /j, e s shown in Fig. [5] and the estimated de- 
viation is only around 6 * 10~ 3 . 

It it astonishing how little difference the polydispersity 
makes. It would clearly be a different situation if we 
had much larger fluctuations in the grain sizes. However, 
the level of polydispersity examined here has basically no 
impact compared to a system of identical beads neither 
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FIG. 7: The ratios of normal stresses are plotted against the 
friction coefficient of the grain-grain contacts. Open symbols 
show a xx /<jyy while closed symbols denote a zz /a yy . The data 
distinguished by circles and squares are obtained in polydis- 
perse and monodisperse systems, respectively. 



on the density nor on the stress tensor. 

The critical state, on the other hand, is very sensitive 
to the inter-particle friction coefficient [i. For example 
when changing fi the critical solid fraction (Fig. [5]) ex- 
plores almost the whole range between the two limits 
^•rlp aud ^rcp that are commonly attributed to the 
random loose and random close packings. 

The fact that the effective friction depends strongly on 
/j, is important in the present scope. This property will 
be exploited later to study the problem of the collective 
rheology. 

We find (Fig. [7]) that the normal stresses in the shear 
plane {a xx and o yy ) remain approximately equal for the 
whole range of \i as we saw this for one particular case in 
section Hi Al Also the normal stress in side direction a zz 
appears to be weaker than the other two. However, the 
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friction \i has a significant influence on this weakening 
effect: While for large [i the weakening is about 14 per- 
cent, it drops down to 3 percent with vanishing friction 
coefficient. 



III. BEYOND THE HOMOGENEOUS CASE 
A. The shear profile 

So far we focused on the homogeneous plane shear and 
examined the properties of the critical state. In this sec- 
tion we turn to the case of an inhomogeneous shear flow 
that is still relatively simple, stationary and quasi static. 
We show that the behavior can not be interpreted based 
on the critical state of the material. 

When our shear cell is filled with one material as in 
the previous section, the flow is indeed homogeneous. We 
find e.g. that the local stress tensor averaged over time is 
the same everywhere in the shear cell. The time averaged 
local velocity of the flow is parallel to the x axis and the 
speed is proportional to the y coordinate. Thus the local 
shear rate is independent of the position. 

The inhomogencity is introduced into the flow by using 
two different materials. The upper half of the shear cell 
(y > Ly/2) is filled with grains of contact friction /i up 
that is chosen to be larger than the value [i\ that is 
used in the lower half (y < L y /2). We recall that there 
is no gravity in the simulation, "up" refers only to the 
orientation of the y axis. Apart from the coefficient of 
the contact friction there is no difference between the two 
materials. 

The main question is how the shear strain is distributed 
throughout the system. The shear cell still remained 
translation invariant in x and z directions and no local 
property is expected to depend on these coordinates. We 
will focus on the y dependence of the measured quanti- 
ties, especially of the local shear strain. For such mea- 
surements the system is divided into thin layers of thick- 
ness Ay parallel to the xz plane. When the value of a 
measured quantity is reported later as the function of y it 
means that it is calculated for each layer, averaged over 
x, z and time over the whole range and over y within the 
given layer. 

To set up and study the above suggested system (or an 
equivalent one) in real experiments might be very diffi- 
cult. However, this shear cell suits well to computer sim- 
ulations and is very instructive regarding the problem of 
collective rheology. As mentioned before, our goal is to 
avoid unnecessary complications and study the simplest 
possible system where the problem can be observed. The 
advantage of our system over traditional shear cells (like 
Couette cell, cylindrical and straight split bottom cells) 
is manifold: (i) the layers that slide next to each other in 
the shear flow are not curved but are straight planes @, 
still stationary flow can be maintained, (ii) Due to the 
high degree of symmetry local quantities arc not multi- 
variate functions. They depend only on one parameter, 
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FIG. 8: (a) Sketch of the shear cell. The material in the 
upper part has larger coefficient of contact friction than that 
in the lower part, (b) and (c) The velocity v(y) and the 
corresponding shear strain is shown, respectively, that are 
suggested by the relation fj, c s(I) obtained for homogeneous 
plane shear, (d) For comparison, the shear strain is drawn 
qualitatively that is observed in the computer simulations. 



the coordinate y. (iii) Due to mechanical equilibrium 
the local shear stress a xv and the normal stress (Tyy are 
bound to be constant (do not depend on y). Thus the 
spatial stress distribution is a priori very simple. 

How does such a system deform? What we could 
naively expect based on the results of the homogeneous 
case is the following. Let us consider the lower mate- 
rial first which is easier to shear. As Oxyl<Jyy does not 
depend on y the apparent effective friction is constant 
throughout the lower material and so is the local inertia 
parameter (Fig. [3]). This gives a constant shear rate "f(y) 
and a linear velocity profile in this region (Fig. [5]). As 
the inertia parameter is kept very small the stress ratio 
&xy/&yy must be close to the plateau value of the effec- 
tive friction of the lower material ^t C rit,io- The same stress 
ratio must be valid also in the upper part of the system. 
However, this stress ratio is not enough to cause shear de- 
formation at any rate in the upper material because here 
the contact friction and therefore also the critical value 
Merit, up is higher than in the lower material (Fig. [6]). Thus 
the overall stress ratio Oxyjayy is smaller than the min- 
imum value (/K C rit,up) that would be needed to maintain 
shear flow. This naive argument leads us to the conclu- 
sion that the upper material does not deform at all. The 
local shear rate j(y) is constant in the lower part and 
zero in the upper one. 

This is in contrast to our numerical measurement 
where different behavior is found (Fig. [5J1) . The observed 
flow does not support the assumption that there is a one 
to one relation between the local effective friction and 
the local inertia parameter. The function /x o ff(-0 that 
has been found in the homogeneous case is clearly not 
valid here. 

Regarding the inhomogeneous case, we present three 
computer simulations A, B and C and the parameters of 
the simulations are as follows. The systems are monodis- 
perse (every grain radius is 1.0), include 5000 grains 
each, and are subjected to a typical global shear strain 
7 g iob = Al shcar/(2ij,) ~ 100. The driving shear rate 'ygiob 
and the vertical pressure a yy , which are kept constant in 
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FIG. 9: The normalized shear strain is shown for systems A 
(circles) and C (squares) as the function of the vertical posi- 
tion y. Both systems have the same position of the interface 
indicated by the double vertical line. The pair of contact fric- 
tions (fj,\ ,fi up ) used below and above the interface is (0.1,0.2) 
for A and (0.1,0.5) for C. The data are fitted by cosine (hy- 
perbolic cosine) functions shown by the solid (dotted) lines. 
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FIG. 10: The normalized shear strain is shown for systems A 
(circles) and B (squares) as the function of the vertical posi- 
tion y. Both systems have the same pair of contact frictions 
(MiojMup) used below and above the interface: (0.1,0.2). They 
differ only in the position of the interface as indicated by the 
double vertical lines. The data are fitted by cosine (hyperbolic 
cosine) functions shown by the solid (dotted) lines. 



time, are chosen to be quasi static: The corresponding 
global inertia parameter / g i b is 5 * 10 -4 in each case. In 
system A and B the pair of contact frictions (ni , A*up) is 
set to (0.1, 0.2) while for system C the contrast is larger: 
(0.1,0.5). In case of A and C the two materials occupy 
half and half of the shear cell as described before. This is 
different for system B where the y-position of the inter- 
face is set to 0.7 L y thus the width of the region of small 
(large) friction is 70 (30) percent of the width L y of the 
shear cell. 

During the shear flow grains of the two materials could, 
in principle, diffuse and mix. This effect is not in the 
present focus, however, in the long run (in the time scale 
of mixing) would significantly alter the shear profile. To 
avoid this interference we exploit the possibility provided 
by the computer simulation and switch off the mixing 
effect as follows. If a grain-grain contact is located below 
or above the interface then its friction coefficient is set 
to fi\ Q and /i up , respectively. Thus contact friction is a 
property of the position and not of the grains in these 
simulations which results in a sharp interface between the 
two materials during the flow. In real experiments mixing 
can not be avoided in this way. We note, however, that 
due to separation of time scales such interfaces can exist 
long enough and mixing does not necessarily become an 
issue [Hill]. 

The data on the observed shear profiles are shown in 
Figs. [9] and [TU] where the local shear strain is plotted 
divided by the global shear strain that was applied to 
the shear cell (7(y)/7 g k>b)' These data are equivalent to 
the normalized shear rate 7(y)/7 g iob and as long as the 
flow is quasi-static the observed curves are independent 
of the driving shear rate 7 g iob- In each system the steady 
state involves shear deformation of both the lower and the 



upper materials. These simulations demonstrate clearly 
that various values of the local shear rate are possible 
for the same local stress. The measured effective friction 
Meff = Cxy/cyy, which can be interpreted as the shear re- 
sistance that the inhomogeneous system collectively de- 
velops against the external driving, turns out to be very 
similar in all cases A, B and C: /i c ff = 0.275 ± 0.003. 
This value is somewhat larger then the estimated critical 
friction for the lower material: jU CT it(A* = 0-1) ~ 0.264 
and smaller than those for the two high-friction materi- 
als: jUcritO" = 0-2) ~ 0.314 and /i cr it(M — 0-5) ~ 0.356. 



B. Tentative description 

How can we interpret the behavior presented in Figs. [5] 
and 1101.'' Why docs, in the first place, the region of ^ U p 
flow at all for shear stresses that are too weak to maintain 
flow in the homogeneous case? Presumably, the material 
of /x up remains solid for such stresses only if there is no 
flow around. However the same material can not retain 
the same static shear stress in the inhomogeneous situa- 
tion because the lower part of the system flows already. 
This flow provides a kind of noise in the force network of 
the upper material. Thus the region that is expected to 
be solid is exposed to permanent agitation which must be 
quite strong near the interface. Under these noisy condi- 
tions the inner fabric might fail from time to time which 
results in a creep motion in the direction of the local 
shear stress no matter how weak this shear stress is. The 
key concept must be in our view that local deformation 
represents also the source of agitation [l(| . Even without 
knowing exactly the nature of the agitation effect it seems 
to be plausible that this mechanism is responsible for the 
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observed quasi static rheology in various systems. The 
agitations arising from the nearby flow was expected to 
rule out static shear stresses also in the experiment [1(3] 
where sand was turned into a zero yield stress fluid. 

Let us set up a tentative model in the spirit of the 
above picture to describe the rheology seen in the previ- 
ous section. We consider an array of narrow parallel lay- 
ers perpendicular to the y axis with their width £ much 
smaller than the system size. The shear strain of the ith 
layer is denoted by 7j. The value a yy , which is a fixed 
constant in the simulation, is taken as the unit of stress 
and we express the shear stress in this unit r = <T xy /a yy . 
We assume about the local shear strain ji that it de- 
pends on two factors, on the local shear stress and on the 
amount of agitation that the ith layer receives. Further- 
more, the deformation of the neighboring layers serves 
as the source of agitation, i.e. the amount of agitation 
received by the layer i is the function of 7»_i and ji+i 
only. We put this in the following form: 



7t = (li-i +7»+i)/(t) , 



(3) 



where the plausible assumption is taken that the caused 
deformation 7, is proportional to the deformation ji—i + 
7f+i that creates it. About the effect of r we assume only 
that it is given by a monotonous increasing function /(r) 
that may be different for different materials. It means 
that for the same amount of agitation the resulting strain 
grows with the shear stress. Eq. [3] is equivalent to the 
discrete equation: 



li-i + 7t+i - 2ji 

e 



2 

e 



1 



2/(r) 



1 



(4) 



In a continuum description (on length scales larger than 
£) this corresponds to the following differential equation 
for 7(2/): 



dy 2 



C(r) 7 , 



(5) 



where C(r) is a material dependent function of the shear 
stress 



C(r) 



2 

e 



1 



2/(r) 



1 



(6) 



The aim of Eq. [5] is to apprehend the shear profile in 
the studied flow which is quasi static and stationary. The 
structural form of the equation seems to be appropriate: 
it does not involve time and time derivatives which is 
a nice feature for a quasi static model, furthermore, if 
7(?/) is a solution then X'j(y) is a solution as well for any 
number A as it should be in steady state. 

Let us consider the homogeneous situation first where 
7(y) is independent of y and t is equal to the critical 
stress ratio T crit . Then the left hand side of Eq. [5] is 
zero and the value of C(r) must vanish for r cr j t . It also 
gives /(T cr it) = 1/2. This has an important consequence 
beyond the homogeneous case, namely, that the sign of 



C is different below and above the critical shear stress. 
Because /(t) is monotonous increasing: 



C{t) > if 
C(t) < if 



T < Tcrit , 
T > Tcrit ■ 



(7) 



Eq. [5] allows quasi static flow also for shear stresses be- 
low and above the critical value T cr it and it makes also 
possible that different parts of the material exhibit dif- 
ferent extent of deformation under the same local stress. 
This is in contrast to the naive picture we deduced alone 
from the homogeneous simulations of section III Al 

The shear resistance that the material sets up against 
the external driving may be smaller or also larger than 
the critical value seen in the homogeneous plane shear. 
Thus not only softening but also hardening is possible 
compared to this reference state. Eq. [5] tells us that 
the deviation from the critical stress results in a curved 
shear profile 7(2/). The curvature is positive (negative) 
for smaller (larger) values of r. 

t is constant in space in our shear cell for both the ho- 
mogeneous and the inhomogencous simulations. There- 
fore C(t) is always a constant throughout one material 
and for such a region the functional form of 7(2/) can be 
easily determined from Eq. [5] 



7io(y) = A\o cos [fci (y - y c )] 



(8) 



for the lower materials in Figs. [91 and [TOl where y c is the 
actual middle position of the given region (where we uti- 
lized the mirror symmetry of the setup with respect to 
the symmetry plane of position y c ). The parameter k\ Q 
is given by C\ (t) = — k 2 Q , where C is indeed negative 
as t > Tcrit and the curvature of 7 is also negative. For 
the upper materials where in each case the overall shear 
stress t is smaller than T cr it of the given material the 
value C up (t) is positive and the corresponding solution 
has positive curvature: 



7u P (y) = -4up cosh [k up (y - y c )} 



(9) 



where fe^ p = C up (t). 

The above model gives nice interpretation of some fea- 
tures of the quasi static behavior. As mentioned already 
both softening and hardening is possible. Hardening is 
interpreted that the ith layer gets less amount of agita- 
tion per unit shear strain than it needs for homogeneous 
deformation. The reduced level of agitation is connected 
to negative curvature of the shear profile 7(2/) and makes 
the material harder to shear (t > T cr it)- Softening com- 
pared to the reference critical state is the opposite case 
which involves enlarged level of agitations, positive cur- 
vature of 7 and smaller shear stress than the critical one. 
The model also does a good job in predicting the func- 
tional form of the shear profile as the functions cosine 
and hyperbolic cosine match the observed data of 7(2/) 
very well (Figs. and [TUl) ■ 

There arc of course many questions left open in our 
quasi one dimensional shear cell which arc important also 
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to understand quasi static flows in general. For exam- 
ple we can not calculate the parameters k\ a and k up as 
we do not know the function C(r). The prediction of 
ratio A\ /A up is also missing here. Therefore all these 
quantities are fit parameters in Figs. l9l and HUl The pre- 
dicted curves arc fitted only inside the given material, 
data points at the boundaries are not included because 
they represent an average shear strain in the vicinity of 
the interface involving both materials. 

In the above picture the shear stress r determines k 
but not the amplitude. In the lower material A\ a de- 
pends only on the amount of agitation that this region 
gets through the boundaries, that is, on the agitation 
that is generated by the other material and transmitted 
through the two interfaces. In turn, A up is determined 
by the transmitted agitation through the interfaces in 
the other direction towards the upper material. There is 
no reason why 7(2/) should be continuous at such mate- 
rial interfaces. At this point we are not able to tell the 
boundary condition for 7(2/) since not enough is known 
about the agitation mechanism. 
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FIG. 11: The local solid fraction is plotted for simulation A 
as the function of the vertical position y. The double vertical 
line indicates the position of the interface between the lower 
and the upper materials. For both materials the horizontal 
dashed lines show the values $ cr it , i.e. the corresponding solid 
fractions in case of a quasi static homogeneous plane shear. 



C. The solid fraction 

Let us consider the inhomogencous shear flow in our 
shear cell and the local state of the material at a given 
y position. How does this state differ from the critical 
state of the same material obtained in the homogeneous 
case? We argued that the two situations are somehow 
different in two quantities, in the amount of the hypo- 
thetical agitation per unit shear strain and in the more 
specific stress ratio G xy j o yy under which the flow occurs. 
Here we would like to emphasize and demonstrate it more 
clearly that the two states are different, even in the inner 
structure of the material. 

This can be well seen if we plot the local solid fraction 
$(?/) measured in the inhomogeneous flow. Fig. [TT] shows 
the data for simulation A, where also the critical solid 
fractions of the two materials are shown for comparison 
that were obtained by homogeneous plane shear. The 
lower material exhibit a solid fraction in this stationary 
and quasi static flow that is smaller than its well defined 
critical value. The upper material does the opposite, it 
gets denser compared to its critical density. We inter- 
pret the phenomenon as follows. The shear deformation 
tends to create pores between the grains while agitations, 
a sort of random shaking, tends to destroy these pores 
and drive the material towards a relatively dense struc- 
ture. Compared to the homogeneous situation the upper 
material is deformed in an environment with increased 
level of agitation therefore the densification works bet- 
ter. For the lower material the reduced level of agitation 
results in looser structure than in the critical state. 

Whether the speculation about the agitation level is 
right or not, the effect is clearly shown that the lower 
material is packed looser than in its critical state while 
the upper material gets denser than its own critical den- 



sity. This effect is so strong that actually the order of 
the densities of the two materials is reversed. This is sur- 
prising because it is expected that spherical grains with 
larger contact friction have lower solid fraction. Here we 
find the opposite within one system in a stationary flow. 



IV. SUMMARY 

This paper is subjected to the quasi static flow of gran- 
ular materials. We argued that quasi static rheology is 
a collective phenomenon that can not be interpreted lo- 
cally as a simple plane shear. In general, the collective 
flow can exhibit such local states of the material that do 
not exist at all in homogeneous plane shear tests. There- 
fore it is difficult to isolate these states and study them 
separately. Even if we restrict the study to quasi static 
and stationary shear flows the physical properties of a 
mesoscopic piece of bulk material can be very different 
from the properties of the well defined critical state (e.g. 
density, shear resistance). The collective rheology can 
not be understood merely based on a constitutive rela- 
tion that connects the local stress to the local strain and 
strain rate. 

We analyzed a rather simple shear cell that was de- 
signed to demonstrate these fundamental problems of col- 
lective granular flows. We used three dimensional com- 
puter simulations and analyzed the spatial distribution 
of the stress and strain. We discussed a physical picture 
that may help to understand the observed flow in which 
the local deformation agitates the surrounding material 
and therefore creates deformation. We set up a tenta- 
tive continuum model that essentially states that the lo- 
cal shear strain 7 is proportional to its Laplacian A7 
and the factor between them is determined by the local 
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stress. The shear profile predicted by the model is in nice search Fund (grant No. PD073172). 
agreement to the observed deformation in our quasi one 
dimensional system. 
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